% Norm

function y = func_norm_res(m, CONSTS, coef_n, p_n)

    c = CONSTS.c;
    a = CONSTS.a;
    
    r_high = 10000;
    num_r_points = 100000;
    r_step = r_high/num_r_points;
    r1 = (0 : r_step : 1)';
    r2 = (1 : r_step : r_high)';
    
    inside  = func_inside(r1, m, CONSTS, coef_n, p_n);
    outside = func_outside(r2, m, CONSTS, coef_n, p_n);
    
%     figure; plot(r1, real(inside), 'b-', r1, imag(inside), 'r-');
%     figure; plot(r2, real(outside), 'b-', r2, imag(outside), 'r-');
    
    inside(isnan(inside)) = 0;
    outside(isnan(outside)) = 0;
    
%     y1 = c*(a^2) * (trapz(r1, inside));
%     y2 = c*(a^2) * (trapz(r2, outside));
    
    y1 = c*(a^2)* (quadgk(@(r)func_inside(r, m, CONSTS, coef_n, p_n), 0.0, 1)); %c*(a^2)*
    y2 = c*(a^2)* (quadgk(@(r)func_outside(r, m, CONSTS, coef_n, p_n), 1, r_high)); %c*(a^2)*
    y = y1 + y2;
    
end

function y = func_inside(r, m, CONSTS, coef_n, p_n)

    [Hrho_in, Hphi_in, Hz_in, Erho_in, Ephi_in, Ez_in] = func_HE_inside(r, m, CONSTS, coef_n, p_n);

    y = (Erho_in .* Hphi_in + Ephi_in .* Hrho_in) .* r;
    
end

function y = func_outside(r, m, CONSTS, coef_n, p_n)

    [Hrho_out, Hphi_out, Hz_out, Erho_out, Ephi_out, Ez_out] = func_HE_outside(r, m, CONSTS, coef_n, p_n);
 
    y = (Erho_out .* Hphi_out + Ephi_out .* Hrho_out) .* r;
    
end

